C.................... CALVN2.FOR;24        12-MAR-1993 21:22:33.86 
C------ This subroutine is the main sequencing control for the vibrational
C------ excited N2 solutions. It was developed in 1993 by splitting up the 
C------ old RSVLPE.FOR. It specifies time and convergence looping and calls
C------ subroutines to set up and solve the coupled continuity and 
C------ momentum equations using Newton's method.
C------ Programmed by Phil Richards in February 1993
C------ Modified by P. Richards in December 2008 to allow for variable
C------ upper boundary altitude that is needed for low L-values. 
C------ Vertical grid cannot exceed the L-shell apex altitude so that 
C------ ion densities can be interpolated to vertical grid
      SUBROUTINE CALVN2(MINJ,MAXJ,S,IDIM,DTV,GLAT,GRD,IC,ISOL
     >   ,MAXIT,IVLPS,HMF2N,HMF2S,SATN,SATF,CHIN,CHIF)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      INTEGER NION,NEQ,NFLAG,VD,LAXIT,JTI
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      REAL AP,DEC,ETRAN,BLON,F107,F107A,SEC,D,GRD,VDUM(VD),RVBV(VD)
      REAL SATN,SATF,CHIN,CHIF
      DIMENSION S(IDIM,1),D(8),VM(7),VK(6)
      COMMON/V91/RVBT(401)
      COMMON/SLV/DELTA(1806),RHS(1806),WORK(50000)
      COMMON/STAW/STR(12,401),RCON(401)
      COMMON/FJS/N(4,401),TI(3,401),FF(20)
      COMMON/NNVV/F(20),NV(7,VD),NVS(7,VD)
      COMMON/ALT/Z(401),DT,DW,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,NION
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GL(401)
      COMMON/NDM/NNO(VD),N2D(VD),O1D(VD),N4S(VD),N2N(VD),ON(VD),CHIV(VD)
      COMMON/ND1/NR(2,VD),ZR(VD),TN(VD),GR(VD),BG(VD),VDT,TE(VD),DH,JTER
      COMMON/ND2/O2N(VD),TIZ(VD),ZNE(VD),HE(VD),N2P(VD),HN(VD)
      COMMON/PXS/VK2(6),PXN2A(VD),PEEX(7,VD),PRN2P(7,VD),JVMAX
      COMMON/RK/VC(401)
      COMMON/LPS/SZA(401),D1,D2,ISPC,I1,I2,IPRIN,JTIMAX,IPP,ISKP
      COMMON/VAGRID/JBV,JVB,JNB
C
       DATA EPSN, RE  , DCR  ,DC, EPCON, NFLAG, IPR, IVN2W, JTI
     >   /   .8 ,6357., 0.0  ,.5, 1.0E-2,  0  ,  6 ,   -1 ,  0/
C------ relative reaction rates N2*(0-5) with O+.. e*-e factors for 0-6
C------ note sixth e*-e factor includes all higher levels
      DATA VK/1,1,38,85,220,270/, VM/7*28./
C
C........ ID=3,4 - output file for northern, southern hemisphere
        ID=3
        IF(MINJ.GT.1) ID=4
        JVMAX=JMAX
        VDT=DTV
        JTER=0
        JTI=JTI+1
        
      !.. Find the upper boundary for vertical drift, making sure it is 
      !.. not greater than the field line apex altitude
      ! Modified in December 2008 to allow for variable upper boundary 
      DO J=1,JBV
        IF(MAXJ.LT.(JMAX+1)/2.AND.ZR(J).LT.Z(MAXJ)) JBP1=J
        IF(MAXJ.GT.(JMAX+1)/2.AND.ZR(J).LT.Z(MINJ)) JBP1=J
      ENDDO
      JBP=JBP1-1
      JBL=JVB           !.. Lower boundary index
      MIT=JBP-JBL+1     !.. # points on the altiude grid to solve
      NEQ=IPR*(MIT-1)   !.. # of equations to solve
C
C-------- Set up terms in the DE that do not change during iteration
      CALL CALDIF(JBP1,1,VM,IC)
C
C----- Generate initial profiles or read from storage STR(IS,J) initial
C----- profiles for N2* are Boltzmann. TINVIB and 1.0E-12 are used to 
C----- ensure that the density doesn't go to zero at low altitudes
      DO 21 J=1,JBP1
         JCON=J
         IF(MINJ.GT.1) JCON=MAXJ+1-J
         DO 20 IS=1,IPR+1
            NV(IS,J)=STR(IS,JCON)
            TINVIB=TN(J)
            !..IF(TINVIB.LT.300) TINVIB=300.0
            IF(ISOL.EQ.0.OR.IS.EQ.IPR+1) NV(IS,J)=10.000+N2N(J)
     >        *DEXP(-(IS-1)*3380.0/TINVIB)/(1.0-DEXP(-3380.0/TINVIB))
 20      CONTINUE
 21    CONTINUE

      !... Find the index of hmF2 on the vertical grid
      DO J=2,JBP1
        IF(MINJ.LE.1) THEN
           IF(HMF2N.LT.(ZR(J+1)+ZR(J))/2.AND.HMF2N.GE.(ZR(J)+ZR(J-1))/2)
     >      JPEAK=J
        ELSE
           IF(HMF2S.LT.(ZR(J+1)+ZR(J))/2.AND.HMF2S.GE.(ZR(J)+ZR(J-1))/2)
     >       JPEAK=J
        ENDIF
      ENDDO

C------ save density at previous time step for calculation of dn/dt
      DO 41 J=1,JBP1
      DO 41 IS=1,IPR
      NVS(IS,J)=NV(IS,J)
 41    CONTINUE
C
C------ Boltzmann N2* boundary condition. TINVIB and E-12 as above
      DO 297 J=1,JBL
         TINVIB=TN(J)
         !..IF(TINVIB.LT.300) TINVIB=300.0
         DO 299 IS=1,IPR
            NV(IS,J)=10.000 + N2N(J)
     >        *DEXP(-(IS-1)*3380.0/TINVIB)/(1.0-DEXP(-3380.0/TINVIB))
 299     CONTINUE
 297  CONTINUE
C
C------ upper boundary density reset to msis value for N2*(v=0) solution
      NV(1,JBP1)=N2N(JBP1)
C
      IF(ISOL.EQ.9) RETURN

C------- MAIN loop - begin iterative solution procedure ends at next row of
C**********************************************************************
      LAXIT=IABS(MAXIT)
      DO 220 JTER=1,LAXIT
C
C----- compute Fij profile and S matrix which BDSLV uses to
C----- solve the linear system ((S))*(DELTA)= (Fij)
      DO 145 J=2,MIT
         JC=J+JBL-1
         CALL VFIJ(ID,MINJ,MAXJ,JC,0,IPR,JBP)
         KR=IPR*(J-2)
         DO 139 IRHS=1,IPR
 139     RHS(KR+IRHS)=F(IRHS)
 145  CONTINUE
C----- IJM is a switch for calling OMATRIX or not ...Cuts down on the amount
C----- of work in building the matrix.
       IJM=0
       IF((JTER/3)*3.EQ.JTER) IJM=1
       IF(JTER.LE.3) IJM=1
       IF(DCR.NE.1) IJM=1
C
       IF(IJM.EQ.1) CALL VMATRX(ID,MINJ,MAXJ,S,RHS,NEQ,IPR,MIT,JBL,JBP)
C
C------- invert the Jacobian matrix S and record error comments
C-------   from the inversion routine BDSLV. IBW = band width
       IBW=2*IPR-1
       CALL BDSLV(NEQ,IBW,S,0,RHS,DELTA,WORK,NFLAG)
       IF(NFLAG.NE.0) THEN
         WRITE(6,193) JTER
 193     FORMAT(1X,'**** Problem in BDSLV, vib N2, JTER=',I3)
         IF(MINJ.LE.1) WRITE(6,*) ' Problem in Northern hemisphere'
         IF(MINJ.GT.1) WRITE(6,*) ' Problem in Southern hemisphere'
         DO 345 J=1,JBP1
           WRITE(6,195) J,ZR(J),(NV(I,J),I=1,6)
 345     CONTINUE
 195     FORMAT(1X,I3,F6.0,1P,22E10.2)
         CALL RUN_ERROR    !.. print ERROR warning in output files
         STOP !3
      ENDIF
C
C----- Update previous densities NV with increment DELTA from BDSLV
C----- Test for convergence when IDIV=0 ... no densities change by more
C----- than specified amount.
 55    IDIV=0
       DCR=1.0
C------- First, test delta to see if it is bigger than NV -> -ve densities
C------- Set DCRP to ensure NV remains positive
      DO 142 I=1,IPR
         ION=2*IPR-I
         DO 142 J=2,MIT
         DCRP=1.0
         DINC=DELTA(IPR*J-ION)
         JC=J+JBL-1
         IF(DINC.LT.0.0) GO TO 142
         IF(ABS(DINC/NV(I,JC)).GT.EPSN) DCRP=DC*ABS(NV(I,JC)/DINC)
         IF(DCRP.LT.DCR) DCR=DCRP
 142  CONTINUE
C
C----- subtract all or part of DELTA from density NV to get new density
C----- for next iteration
 144  DO 42 J=2,MIT
         JC=J+JBL-1
         DO 42 I=1,IPR
         ION=2*IPR-I
         DINC=DELTA(IPR*J-ION)
         NV(I,JC)=NV(I,JC)-DINC*DCR
         IF(ABS(DINC/NV(I,JC)).GT.EPCON)  IDIV=IDIV+1
 42   CONTINUE
C
C------ If all densities have a small enough change, solution achieved
 47   IF(IDIV.EQ.0.AND.JTER.GT.17) WRITE(6,113) JTER
 113  FORMAT(/1X,'ITER=',I5,' ** Difficulty converging in N2* soln'//)
      IF(IDIV.EQ.0)  GO TO 222
C
 220  CONTINUE

      WRITE(6,*)
      WRITE(6,*) ' ****** non-convergence in vibrational N2'
      CALL RUN_ERROR    !.. print ERROR warning in output files
      STOP

C**********************************************************************
C
 222  CONTINUE

C----- calculate new rate constant VKEN, NTOT, VKE are used in
C----- calculation. see Newton et al p3813 eqn 27
      DO 88 J=1,JBP1
         JCON=J
         IF(MINJ.GT.1) JCON=MAXJ+1-J
         VKE=0.0
         NTOT=0.0
C------ calculate reduction in N2 cooling rate of thermal electrons
         CALL STUBV(2,TE(J),TN(J),ZNE(J),N2N(J),RBOLT,NV,J)
         CALL STUBV(1,TE(J),TN(J),ZNE(J),N2N(J),RVIB,NV,J)
         IF(RBOLT.NE.0) RVBV(J)=RVIB/RBOLT
         DO 223 I=1,IPR
            D(I)=SNGL(NV(I,J))
            STR(I,JCON)=NV(I,J)
            NTOT=NTOT+NV(I,J)
C------ calculate Boltzmann distribution for comparison
C...      IF(JTIMAX.EQ.1) NVT(I)=NV(1,J)
C...     >  *DEXP(-(I-1)*3380.0/TN(J))/(1.0-DEXP(-3380.0/TN(J)))
 223        VKE=VKE+VK(I)*NV(I,J)
         VKEN=VKE/NTOT
         VDUM(J)=SNGL(VKEN/VK(1))
         STR(7,JCON)=VKEN/VK(1)
         IF(J.GT.0) GO TO 88
            JP1=J-1
            IF(J.EQ.1) JP1=4
            TV=TN(J)
            IF(ZR(J).GE.ZR(JVB)) CALL TVIB(J,TN(J),TV,D)
            IF(J.EQ.1) WRITE(6,117)
 117        FORMAT(/'  ALT   TV   TN    VKEN',3X,'(N2:',5X,'V=0'
     >        ,5X,'V=1',5X,'V=2',5X,'V=3',5X,'V=4',5X,'V=5')
            IF(ZR(J).LT.100)
     >      WRITE(6,116) ZR(J),TV,TN(J),VKEN,N2N(J),(NV(IS,J),IS=1,IPR)
 88   CONTINUE

      !.. At low lats, put values between the last vertical and field line grid pts
      DO J=1,JMAX
         IF(Z(JMAX/2+1).LT.600.AND.Z(J).GE.ZR(JBP)) THEN
         RCON(J)=RCON(J-1)
c         WRITE(217,'(I4,9F8.1)') J,ZR(JBP),ZR(JBP1),Z(J),RCON(J)
c     >    ,RCON(J-1),VC(J-1),VC(J-4),VC(J-7),VC(J-9)
         ENDIF
      ENDDO

       !... Write N2* density and temp at hmF2
 119   FORMAT(/3X,' ..... Vibrational N2 temperature and populations'
     >   ,1X,'at hmF2 ....'
     >   ,/6X,'UT    LT  CHI  ALT   TV    TN  kN2v  MSISN2'
     >   ,4X,'v_0    v_1     v_2     v_3     v_4      v_5')
       IF(IVN2W.EQ.1) THEN
         IF(JTI.EQ.1) WRITE(40,119)
         IF(JTI.EQ.1) WRITE(41,119)
         IF(JPEAK.GT.1) THEN
           J=JPEAK
           DO I=1,IPR
             D(I)=SNGL(NV(I,J))
          ENDDO
          TV=TN(J)
          IF(ZR(J).GE.ZR(JVB)) CALL TVIB(J,TN(J),TV,D)
          IF(MINJ.LE.1) THEN
            WRITE(40,118) SEC/3600.,SATN,NINT(CHIN*57.2958),NINT(ZR(J))
     >        ,NINT(TV),NINT(TN(J)),VC(J),N2N(J),(NV(IS,J),IS=1,IPR)
          ELSE
            WRITE(41,118) SEC/3600.,SATF,NINT(CHIF*57.2958),NINT(ZR(J))
     >        ,NINT(TV),NINT(TN(J)),VC(J),N2N(J),(NV(IS,J),IS=1,IPR)
          ENDIF
         ENDIF
       ENDIF 
      
C------interpolate rate constant factor 
      DO 955 J=1,JMAX/2
        JCON=J
        IF(MINJ.GT.1) JCON=MAXJ+1-J
        RVBT(JCON)=1.0
        RCON(JCON)=1.0
        VC(JCON)=1.0
        IF(J.GT.VD.OR.Z(JCON).GT.ZR(JBP).OR.Z(JCON).LT.ZR(1)) GO TO 955
          VC(JCON)=FTVC(1,JBP1,ZR,VDUM,Z(JCON))
          RCON(JCON)=VC(JCON)
          RVBT(JCON)=FTVC(1,JBP1,ZR,RVBV,Z(JCON))
          IF(TI(3,JCON).LT.500.OR.TI(3,JCON).GT.10000) RVBT(JCON)=1.0
          IF(J.LT.0)
     >     WRITE(6,205) JCON,Z(JCON),RVBT(JCON),RVBV(JCON),RCON(JCON)
     >      ,STR(1,JCON),STR(2,JCON),STR(3,JCON)
 955  CONTINUE
 965  CONTINUE

c      IF(MAXJ.LT.(JMAX+1)/2.) THEN
c       JPEAK=JPEAK
c       CALL CHEMO_VN2(6,IPR,JPEAK,ZR(JPEAK),N2N(JPEAK),NV,VKNEW)
c       WRITE(81,'(5F9.2,1P,9E10.2)') SEC/3600.,10.0+SEC/3600.,SATF
c     >  ,CHIF*57.3,ZR(JPEAK),VDUM(JPEAK),VKNEW,VDUM(JPEAK)/VKNEW
c       JPEAK=JPEAK
c      ENDIF
      !--- write individual cpts of cty eqn in VFIJ, or PRLN if IV<50.
      IF(IVLPS.GT.0.OR.IPP.EQ.0) RETURN
      DO 93 IV=2,5
      DO 93 J=1,JBP,IPP
         IF(IV.NE.4.AND.IV.NE.5) GO TO 93
         IF(J.GT.1) CALL VFIJ(ID,MINJ,MAXJ,J,IV,IPR,JBP)
 93   CONTINUE
       RETURN
C
 116  FORMAT(F5.0,2F6.0,F5.2,1P,22E8.1)
 118  FORMAT(F8.2,F6.2,I4,I5,2I6,F5.2,1P,22E8.1)
 205  FORMAT(I5,F7.0,1P,22E9.2)
      END
C::::::::::::::::::::::: VFIJ ::::::::::::::::::::::::::::::::
        SUBROUTINE VFIJ(ID,MINJ,MAXJ,J,ILJ,IPR,JBV)
C------ this program solves the diffusion-continuity eqns for the first ipr
C------ states of vibrational n2. for the theory, see newton et al  1974
C------ jgr p 3817. for constant step lengths put all bg(j)=1 and set
C------ dh=the desired step length in cm. this program sets up the error
C------ functions which are solved by the newton method.
C------ these programs were written by p richards in september 1979
      IMPLICIT DOUBLE PRECISION(A-H,K-L,N-Z)
      INTEGER VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      DIMENSION B(22),FYUBV(6)
      COMMON/DD/D(7,VD),E(7,VD),B1(7,VD),B2(7,VD),RDZ(VD),EDDY
      COMMON/NNVV/F(20),NV(7,VD),NVS(7,VD)
      COMMON/ALT/Z(401),DT,DW,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,IION
      COMMON/NDM/NNO(VD),N2D(VD),O1D(VD),N4S(VD),N2N(VD),ON(VD),CHIV(VD)
      COMMON/ND1/NR(2,VD),ZR(VD),TN(VD),GR(VD),BG(VD),VDT,TE(VD),DH,JTER
      COMMON/ND2/O2N(VD),TIZ(VD),ZNE(VD),HE(VD),N2P(VD),HN(VD)
      COMMON/PXS/VK(6),PXN2A(VD),PEEX(7,VD),PRN2P(7,VD),JVMAX
      DATA FYUBV,IR69/6*0.0,0/
C
      STS=1.0
      IF(ITF.GE.3) STS=0.0
C
      DO 20 IS=1,IPR
C------upper boundary condition.......
      IF(J.EQ.JBV) NV(IS,J+1)=-E(1,J)*NV(IS,J)/RDZ(J)
     >     +NV(IS,J-1)-FYUBV(IS)/RDZ(J)/D(1,J)
C------coefficients of the differential equation......
      A1=(D(1,J)+EDDY)*(BG(J)/DH)**2*(NV(IS,J+1)-2*NV(IS,J)+NV(IS,J-1))
      A2=B1(1,J)*(NV(IS,J+1)-NV(IS,J-1))*RDZ(J)
      A3=B2(1,J)*NV(IS,J)
      CALL PRLO(MINJ,MAXJ,ILJ,IS,J,IPR,PROD,LOSS,B)
C
C----- f(is)= is the sum of the terms in the continuity equation which
C----- we are trying to minimize
      F(IS)=(A1+A2+A3)-(NV(IS,J)-NVS(IS,J))*STS/VDT+(PROD-LOSS)
C
C------ write various cpts of cty eqn........
      IF(ILJ.LE.1) GO TO 20
      A4=(NV(IS,J+1)-NV(IS,J-1))*RDZ(J)
      A5=E(1,J)*NV(IS,J)
      FLUX=-D(1,J)*(A4+A5)
      IF(ILJ.EQ.5.AND.IS.EQ.1) WRITE(ID,101)
 101  FORMAT(/2X,'J',4X,'ALT',6X,'FLUX',6X,'DIVF',6X,'TPROD',6X,'TLOSS',
     >  6X,'DN/DT',8X,'A1',9X,'A2',9X,'A3',9X,'A4',9X,'A5',9X,'F')
      DIVF=A1+A2+A3
      DNDT=(NV(IS,J)-NVS(IS,J))*STS/VDT
      IF(ILJ.EQ.5) WRITE(ID,100) J,ZR(J),FLUX,DIVF,PROD,LOSS,DNDT
     >  ,A1,A2,A3,A4,A5,F(IS)
 20   CONTINUE
 100  FORMAT(I4,F6.0,1P,19E11.2)
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      BLOCK DATA BLK6
      IMPLICIT REAL(A-H,L,N-Z)
      COMMON/PEE/P(4,6,6),G(6),HCOK
      COMMON/DPEE/DP(4,6,6)
C------ the following data(g @ pc) are related to thermal electron excitation
C------ of vib n2. the g's are used to calculate the dexcitation.
C------ the p parameters are given by newton et al p3817 for electron exCit.
      DATA G/0.0,2345.2,4632.6,6905.5,9150.,11703.6/ ,HCOK/1.439/
      DATA P/144*0.0/
      DATA (DP(IP,2,1),IP=1,4)/-1.8374E1,-2.2185E4, 1.6843E7,-4.5588E9/
      DATA (DP(IP,3,1),IP=1,4)/-1.8871E1,-1.8510E4,-1.3173E6, 2.3103E9/
      DATA (DP(IP,6,1),IP=1,4)/-2.0258E1,-1.6288E4,-2.5090E6, 7.7905E8/
      DATA (DP(IP,6,2),IP=1,4)/-2.0258E1,-1.6288E4,-2.5090E6, 7.7905E8/
      DATA (DP(IP,2,3),IP=1,4)/-1.8748E1,-1.2227E4, 1.1025E7,-2.9931E9/
      DATA (DP(IP,6,3),IP=1,4)/-1.9748E1,-1.3149E4,-2.7395E6, 9.1671E8/
      DATA (DP(IP,1,4),IP=1,4)/-1.9504E1,-5.7050E3,-6.3620E6, 1.9172E9/
      DATA (DP(IP,2,4),IP=1,4)/-1.9724E1,-4.6525E3,-6.5852E6, 1.9314E9/
      DATA (DP(IP,3,4),IP=1,4)/-1.9397E1,-8.0303E3, 5.6072E6,-1.2575E9/
      DATA (DP(IP,6,4),IP=1,4)/-1.9807E1,-8.8002E3,-2.5262E6, 9.7401E8/
      DATA (DP(IP,1,5),IP=1,4)/-1.9641E1,-3.9130E3,-5.7353E6, 1.8426E9/
      DATA (DP(IP,2,5),IP=1,4)/-2.0194E1,-3.5886E3,-4.8620E6, 1.6888E9/
      DATA (DP(IP,3,5),IP=1,4)/-1.9997E1,-4.4211E3,-3.1949E6, 1.2434E9/
      DATA (DP(IP,4,5),IP=1,4)/-1.9438E1,-4.8694E3, 2.0787E6,-1.9350E8/
      DATA (DP(IP,6,5),IP=1,4)/-1.9832E1,-4.7005E3,-4.2176E5, 4.7237E8/
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE INITP(IFRST)
      IMPLICIT REAL(A-H,L,N-Z)
      COMMON/PEE/P(4,6,6),G(6),HCOK
      COMMON/DPEE/DP(4,6,6)
      IFRST=1
      DO 10 IP=1,4
        P(IP,2,1)=DP(IP,2,1)
        P(IP,3,1)=DP(IP,3,1)
        P(IP,6,1)=DP(IP,6,1)
        P(IP,6,2)=DP(IP,6,2)
        P(IP,2,3)=DP(IP,2,3)
        P(IP,6,3)=DP(IP,6,3)
        P(IP,1,4)=DP(IP,1,4)
        P(IP,2,4)=DP(IP,2,4)
        P(IP,3,4)=DP(IP,3,4)
        P(IP,6,4)=DP(IP,6,4)
        P(IP,1,5)=DP(IP,1,5)
        P(IP,2,5)=DP(IP,2,5)
        P(IP,3,5)=DP(IP,3,5)
        P(IP,4,5)=DP(IP,4,5)
        P(IP,6,5)=DP(IP,6,5)
 10     CONTINUE
        RETURN
        END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE PRLO(MINJ,MAXJ,ILJ,I,JP,IPR,PROD,LOSS,B)
C--------this subroutine evaluates the production and loss for the variou
C--------states of vib n2. see newton et al jgr 1974 p3817.
      IMPLICIT REAL(A-H,L,N-Z)
      INTEGER VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      DOUBLE PRECISION PROD,LOSS,B,STR,RCON,F,NV,NVS,NNO,N2D,O1D,N4S,
     > ON,NR,ZR,TN,GR,BG,VDT,TE,VK,PXN2A,PEEX,PRN2P,O2N,TIZ,ZNE,HE,N2P
     > ,HN,RTS,N2N,CHIV,DH
      DOUBLE PRECISION TE2,TE3,AIJ,AJI,AIJN,AJIN,VCON,SQRTN,NU,P1001,
     > P1NU,EXPTN,C01,C01ON,TE1
      DIMENSION B(22),RTS(199)
      COMMON/STAW/STR(12,401),RCON(401)
      COMMON/NNVV/F(20),NV(7,VD),NVS(7,VD)
      COMMON/NDM/NNO(VD),N2D(VD),O1D(VD),N4S(VD),N2N(VD),ON(VD),
     > CHIV(VD)
      COMMON/ND1/NR(2,VD),ZR(VD),TN(VD),GR(VD),BG(VD),VDT,TE(VD),DH,
     > JTER
      COMMON/PXS/VK(6),PXN2A(VD),PEEX(7,VD),PRN2P(7,VD),JVMAX
      COMMON/ND2/O2N(VD),TIZ(VD),ZNE(VD),HE(VD),N2P(VD),HN(VD)
      COMMON/PEE/P(4,6,6),G(6),HCOK
      DATA IFRST/0/
C------ the following data(g @ pc) are related to thermal electron excitation
C------ of vib n2. the g's are used to calculate the dexcitation.
C------ the p parameters are given by newton et al p3817 for electron exCit.
C
C----- the b's are used to sum the v-v exchange between diff states. however
C----- b(1) @ b(2) are used to sum elec excit. nu=coll freq eqtn(18). p1001=
C----- dimensionless exchange transition probability for n2 trans 1<--->0
C----- c01= tras-vib energy exchange prob. for colls between o and n2. however
C----- the results from mcneal et al jgr 1974 are used for reaction rate
      IF(IFRST.EQ.0)CALL INITP(IFRST)
      DO  5 IS=1,22
  5   B(IS)=0.0
      CALL RATS(JP,TE(JP),TIZ(JP),TN(JP),RTS)
      SQRTN=DSQRT(TN(JP))
      NU=1.482094E-11*SQRTN
      P1001=2.600322E-6*TN(JP)*EXP(91.5/TN(JP))
      P1NU=P1001*NU
      EXPTN=EXP(3380/TN(JP))
      C01=TN(JP)**2.87*2.4685E-22/EXPTN
      C01ON=C01*ON(JP)
      TE1=1.0/TE(JP)
      IF(TE(JP).LT.400) TE1=1.0/400.0
      TE2=TE1*TE1
      TE3=TE2*TE1
C
C++++++ compts of prod and loss eqtns 16 @ 17 of newton. note that indices
C++++++ in coefficients of production and loss terms
C+++++++ appear reduced by 1, because the fortran indices start from 1, not 0
      DO 40 J=1,IPR
      IF(I.EQ.J) GO TO 50
      IF(P(1,I,J).NE.0.0)
     >  AIJ=EXP(P(1,I,J)+P(2,I,J)*TE1+P(3,I,J)*TE2+P(4,I,J)*TE3)
      IF(P(1,J,I).NE.0.0)
     >  AJI=EXP(P(1,J,I)+P(2,J,I)*TE1+P(3,J,I)*TE2+P(4,J,I)*TE3)
      AIJN=AIJ
      AJIN=AJI
      IF(P(1,I,J).EQ.0) AIJN=AJI*EXP(+(G(J)-G(I))*HCOK*TE1)
      IF(P(1,J,I).EQ.0) AJIN=AIJ*EXP(+(G(I)-G(J))*HCOK*TE1)
C----- zne(jp) contains ne from sub prog prln
      IF(ILJ.NE.4) B(1)=B(1)+AIJN*NV(J,JP)*(ZNE(JP))
      IF(ILJ.NE.4) B(2)=B(2)+AJIN*(ZNE(JP))*NV(I,JP)
      IF(ILJ.EQ.4.AND.J.LT.I) B(1)=B(1)+AIJN*NV(J,JP)*(ZNE(JP))
      IF(ILJ.EQ.4.AND.J.LT.I) B(2)=B(2)+AJIN*(ZNE(JP))*NV(I,JP)
 50   IF(I.EQ.1) GO TO 60
      IF(J+1.GT.IPR) GO TO 60
C------ if i-j=1 quanta are exchanged but no net change in any density
      IF(I-J.EQ.1) GO TO 60
      B(5)=B(5)+P1NU*(I-1)*J*NV(I-1,JP)*NV(J+1,JP)
      B(3)=B(3)+P1NU*(I-1)*J*NV(J,JP)*NV(I,JP)
 60   IF(J.EQ.1) GO TO 40
      IF(I+1.GT.IPR) GO TO 40
C------ if j-i=1 quanta are exchanged but no net change in any density
      IF(J-I.EQ.1) GO TO 40
      B(4)=B(4)+P1NU*(J-1)*I*NV(J,JP)*NV(I,JP)
      B(6)=B(6)+P1NU*(J-1)*I*NV(I+1,JP)*NV(J-1,JP)
 40   CONTINUE
C
C
C******* b(7)=prod due to t-v between o and n2
      IF(I.NE.1) B(7)=C01ON*((I-1)*NV(I-1,JP)+I*EXPTN*NV(I+1,JP))
      IF(I.EQ.1) B(7)=C01ON*I*EXPTN*NV(I+1,JP)
C******* b(8)=loss due to quenching of n2 by o
      IF(I.LE.9) B(8)=C01ON*((I-1)*EXPTN+I)*NV(I,JP)
C******* b(9)=prod thru no+n-->n2*+o & o(1d)+n2--->. & n2(a)
C------ for pno @ po1d see waite et al pss 1979 p901.
      IF(I.EQ.5) B(9)=RTS(9)*N4S(JP)*NNO(JP)
      IF(I.EQ.3) B(9)=RTS(33)*O1D(JP)*NV(1,JP)
C------ production due to the electronic excitation of n2(a)
      IF(I.EQ.1) B(15)=1.3*PXN2A(JP)*NV(1,JP)
      IF(I.EQ.IPR) B(13)=1.3*PXN2A(JP)*NV(1,JP)
C------ quenching of n2* by o2 .......
C...      if(i.ne.1) b(10)=1.0e-13*o2n(jp)*nv(i,jp)
C...      b(11)=1.0e-13*o2n(jp)*nv(i+1,jp)
C------ loss by o+ + n2* --> o + n2+
      JCON=J
      IF(MINJ.GT.1) JCON=MAXJ+1-J
      VCON=STR(7,JCON)
      B(12)=VCON*RTS(3)*NV(I,JP)*NR(1,JP)
C------ production from N2+(v) from subr. VIBN2P. Not being used at
C------ present
C....      IF(I.NE.1) B(13?)=PRN2P(I,JP)
C------ direct excitation by photoelectrons
      IF(I.NE.1) B(14)=PEEX(I,JP)*NV(1,JP)
      IF(I.EQ.1) B(16)=PEEX(I,JP)*NV(1,JP)
C
C------ sum the cmpts of prod and loss
      PROD=B(1)+B(5)+B(6)+B(7)+B(9)+B(11)+B(13)+B(14)
      LOSS=B(2)+B(3)+B(4)+B(8)+B(10)+B(12)+B(15)+B(16)
      IF(I.EQ.1) PROD=0.0
      IF(I.EQ.1) LOSS=0.0
C------ adjust 'loss' for boundary condition calculation
      IF(ILJ.EQ.-9) THEN
        LOSS=LOSS/NV(I,JP)
      ENDIF
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
        SUBROUTINE PRVFIJ(ID,MINJ,MAXJ,J,ILJ,IPR,JBV,IC,VCON)
      IMPLICIT DOUBLE PRECISION(A-H,K-L,N-Z)
      INTEGER VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      DIMENSION B(22),VCON(VD),VPROD(6),VLOSS(6)
      COMMON/NNVV/F(20),NV(7,VD),NVS(7,VD)
      COMMON/STAW/STR(12,401),RCON(401)
      COMMON/NDM/NNO(VD),N2D(VD),O1D(VD),N4S(VD),N2N(VD),ON(VD),
     >  CHIV(VD)
      COMMON/ND1/NR(2,VD),ZR(VD),TN(VD),GR(VD),BG(VD),VDT,TE(VD),DH,
     >  JTER
      COMMON/ND2/O2N(VD),TIZ(VD),ZNE(VD),HE(VD),N2P(VD),HN(VD)
      COMMON/PXS/VK(6),PXN2A(VD),PEEX(7,VD),PRN2P(7,VD),JVMAX
C
      WRITE(ID,102)
 102  FORMAT(/4X,'ALT',5X,'E*',5X,'E-PROD',3X,'E-LOSS',3X,'B3-LOSS',2X
     > ,'B4-LOSS',2X,'B5-PROD',2X,'B6-PROD',3X,'O-PROD',3X,'O-LOSS',1X
     > ,'N+NO&O1D',1X,'O2-LOSS',2X,'O2-PROD',2X,'O+ +N2',5X,'N2A')
C
      JCON=J
      IF(MINJ.GT.1) JCON=MAXJ+1-J
      DO 13 IS=1,4
         NV(IS,J)=STR(IS+7,JCON)
 13   CONTINUE
C
      N2D(J)=STR(8,JCON)
      N4S(J)=STR(9,JCON)
      NNO(J)=STR(10,JCON)
      O1D(J)=STR(11,JCON)
      IF(O1D(J).LE.0) CALL PRLN(0,0,J,0,3,VPROD,VLOSS,VCON,0)
C
      DO 15 IS=1,IPR
         NV(IS,J)=STR(IS,JCON)
 15   CONTINUE
C
      DO 20 IS=1,IPR
      CALL PRLO(MINJ,MAXJ,ILJ,IS,J,IPR,PROD,LOSS,B)
C------ transfer b(15) to b(9), and b(16) to b(14) for print
      IF(IS.EQ.1) B(14)=B(16)
      WRITE(ID,99) ZR(J),B(14),(B(I),I=1,13)
 20   CONTINUE
 99   FORMAT(F6.1,1P,14E9.2)
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE VMATRX(ID,MINJ,MAXJ,S,RHS,NEQ,IPR,MIT,JBL,JBV)
C---- this program sets up the jacobian matrix for the band solver
      IMPLICIT DOUBLE PRECISION(A-H,N-Z)
      INTEGER NEQ,VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      DIMENSION RHS(1806),S(NEQ,23)
      COMMON/NNVV/F(20),NV(7,VD),NVS(7,VD)
      IBW=4*IPR-1
C
      DO 180 JZS=1,NEQ
      DO 190 KZS=1,IBW
      S(JZS,KZS)=0.0
  190 CONTINUE
  180 CONTINUE
C
      DO 80 JF=2,MIT
      J1=MAX0(2,JF-1)
      J2=MIN0(JF+1,MIT)
      DO 90 JV=J1,J2
      DO 130 IV=1,IPR
C
      L=IPR*(JF-2)
      IF(JF.LE.3)M=IPR*(JV-2)+IV
      IF(JF.GT.3)M=IPR*(JV-JF)+IV+2*IPR
      JVC=JBL+JV-1
      JFC=JBL+JF-1
C
      H=1.0E-2*NV(IV,JVC)
      NV(IV,JVC)=NV(IV,JVC)+H
      CALL VFIJ(ID,MINJ,MAXJ,JFC,1,IPR,JBV)
      NV(IV,JVC)=NV(IV,JVC)-H
      DO 75 IS=1,IPR
      IF(JF.LE.3) S(L+IS,M)=(F(IS)-RHS(L+IS))/H
      IF(JF.GT.3) S(L+IS,M-IS)=(F(IS)-RHS(L+IS))/H
 75    CONTINUE
C
  130 CONTINUE
   90 CONTINUE
   80 CONTINUE
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE STUBV(IS,TE,TN,NE,N2,RLVN2,RNV,JP)
      INTEGER VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      DOUBLE PRECISION TE,TN,NE,N2,RLVN2,RNV(7,VD),B
      DOUBLE PRECISION TDIF,EF,GE,TE1,TE2,TE3,AIJ,AJI,AIJN,AJIN,NV
      DIMENSION NV(7),B(2)
      COMMON/PEE/P(4,6,6),G(6),HCOK
C+++ stubbe and varnum  n2 vib loss rates from schunk & nagy p364  +++++
        RLVN2=0.0
      IF(TE.LT.500.OR.TE.GT.10000) RETURN
        TDIF=TE-TN
        EF=1.06E+4+7.51E+3*TANH(1.1E-3*(TE-1800.0))
        GE=3300+1.233*(TE-1000.0)-2.056E-4*(TE-1000.0)*(TE-4000.0)
        RLVN2=-2.99E-12*NE*N2*EXP(EF*(TE-2000)/
     >   (2000*TE))*(EXP(-GE*TDIF/(TE*TN))-1)
C...        if(is.eq.1) rlvn2=3.33*rlvn2
C
C------ now the newton et al loss rates. the p's are given on page 3817
C------ using boltzmann distribution of n2* with temp tn
      DO 35 I=1,6
      IF(IS.EQ.2) NV(I)=N2*EXP(-(I-1)*3380/TN)/(1-EXP(-3380/TN))
      IF(IS.NE.2) NV(I)=RNV(I,JP)
 35   CONTINUE
C
      B(1)=0.0
      B(2)=0.0
      TE1=1/TE
      TE2=TE1*TE1
      TE3=TE2*TE1
C
      DO 40 I=1,6
      DO 40 J=1,6
      IF(I.EQ.J) GO TO 40
      IF(P(1,I,J).NE.0.0)
     >  AIJ=EXP(P(1,I,J)+P(2,I,J)*TE1+P(3,I,J)*TE2+P(4,I,J)*TE3)
      IF(P(1,J,I).NE.0.0)
     >  AJI=EXP(P(1,J,I)+P(2,J,I)*TE1+P(3,J,I)*TE2+P(4,J,I)*TE3)
      AIJN=AIJ
      AJIN=AJI
      IF(P(1,I,J).EQ.0) AIJN=AJI*EXP(+(G(J)-G(I))*HCOK*TE1)
      IF(P(1,J,I).EQ.0) AJIN=AIJ*EXP(+(G(I)-G(J))*HCOK*TE1)
C
      IF(J.LT.I) B(1)=B(1)+AIJN*NV(J)*NE*(I-J)
      IF(J.LT.I) B(2)=B(2)+AJIN*NE*NV(I)*(I-J)
 40   CONTINUE
      RLVN2=B(1)-B(2)
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE TVIB(J,TN,TV,D)
C------ to find the vibrational temperature of the equivalent Boltzmann
C------ N2* densities. Add up the energy in the
C------ calculated vibrational distribution then calculate the vibration
C------ temperature of the Boltzmann distribution that has the same ener
C------ A Newton procedure is employed to make E(Bolt)=E(cal). TN= neutr
C------ temperature for first quess, TV=calculated vibrational temperatu
C------ D(8) contains the calculated N2+(v=0->5). Written by P. Richards
      DOUBLE PRECISION TN,TV,H,DEX,SUM,FEX,FUN
      DIMENSION D(8)
C------ total energy in calculated distribution = D1 E1 + D2 E2 +....
C------ E2=2*E1, E3=3*E1 etc. E1 also appears in Bolt. and cancels out
      SUM=0.0D0
      SUM=D(2)+2*D(3)+3*D(4)+4*D(5)+5*D(6)
C------ first guess=Tn, H= change in TV used to calculate derivative
      TV=TN
      JITER=0
      H=0.0001*TV
C------ find function value for TV
 20   CALL TFUN(TV,D(1),SUM,FUN)
      FEX=FUN
      TV=TV+H
C------ find function value for TV=TV+H for derivative DEX
      CALL TFUN(TV,D(1),SUM,FUN)
      DEX=(FUN-FEX)/H
      IF(DEX.EQ.0) GO TO 30
C------ increment TV
      TV=TV-H-FEX/DEX
      JITER=JITER+1
      IF(JITER.GT.22) GO TO 30
C------ if DEX small enough, stop.
      IF(DABS(FEX/(DEX*TV)).GT.1.0E-3) GO TO 20
      RETURN
C------ something has gone wrong !!
 30   CONTINUE
      TV=0.0D0
      RETURN
      END
C:::::::::::::::::::::::::::::::: TFUN :::::::::::::::::::::::::::::::::
      SUBROUTINE TFUN(TV,D1,SUM,FUN)
C------ Calculate the energy in the Boltzmann distribution temp=TV.
C------ SUM=energy in calculated distribution, D1=density of v=0
C------ FUN= E(Bolt)-E(Cal)
      DOUBLE PRECISION X,TV,X2,X3,X4,X5,BSUM,SUM,FUN
      X=DEXP(-3380./TV)
      X2=X*X
      X3=X2*X
      X4=X3*X
      X5=X4*X
      BSUM=D1*(X+2*X2+3*X3+4*X4+5*X5)/(1.0-X)
      FUN=BSUM-SUM
      RETURN
      END
C::::::::::::::::::::::::::::: FTVC ::::::::::::::::::::::::::::::::::::::
      FUNCTION FTVC(I,JMAX0,ZR,T,Z)
      INTEGER VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      DOUBLE PRECISION ZR,Z,FTVC      ! f90 added
      DIMENSION ZR(VD),T(VD)
      
      DO K=1,JMAX0
        IF(Z.GE.ZR(K).AND.Z.LT.ZR(K+1)) GO TO 2
      ENDDO

 2     FTVC=((Z-ZR(K))/(ZR(K+1)-ZR(K))*(T(K+1)-T(K))+T(K))

      RETURN
      END
C::::::::::::::::::::::::::::::::: BLOCK DATA BLK5 ::::::::::::::::::::::
      BLOCK DATA BLK5
      IMPLICIT DOUBLE PRECISION (A-H,L,N-Z)
      INTEGER VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      COMMON/PXS/VK(6),PXN2A(VD),PEEX(7,VD),PRN2P(7,VD),JVMAX
      DATA VK/1.0,1.0,38.0,85.0,220.0,270.0/
      END
